namespace spectrum {

      using namespace msl;


      typedef enum(
            Base,
            Contig,
            Cyclic,
            Free ) Dist_t;

      class DistributionBase {
      
            protected :
            
                  // Distribution type
                  Dist_t type;
                  // Number of processes
                  int np;
                  // Total rows of DistributedStructuredMatrix
                  int n;
                  // Total columns
                  int m;
                  
                  // Number of submatrix blocks per row
                  int submatricesPerRow;
                  // Number of submatrix blocks per col
                  int submatricesPerCol;
                  // Total submatrices = submatricesPerRow * submatricesPerCol
                  int totalSubmatrices;
                  
                  // Processor ids, ordered by submatrix id
                  std::vector< int > procPerSubmatrix;
                  // Submatrix ids, ordered by processor id and then by submatrix id
                  std::vector< std::vector< int > > submatricesPerProc;
                  // If matrix is contiguous or cyclic, = procPerSubmatrix[0] - 0
                  int offset;

                  // Boundaries between row blocks, in increasing order by row
                  std::vector< int > rowBreaks;
                  // Boundaries between column blocks, in increasing order by col
                  std::vector< int > colBreaks;
                  
                  // Submatrix row extents, ordered by submatrix id
                  std::vector< int > rowSizes;
                  // Submatrix column extents, ordered by submatrix id
                  std::vector< int > colSizes;
                  // Submatrix row boundaries, ordered by submatrix id
                  std::vector< std::pair< int, int > > rowRanges;
                  // Submatrix column boundaries, ordered by submatrix id
                  std::vector< std::pair< int, int > > colRanges;
                  
                  // Submatrix block numbers, ordered by submatrix id
                  std::vector< std::pair< int, int > > blockNumbers;
                  // Submatrix ids stored by block
                  std::vector< std::vector< int > > submatrixPerBlock;
                  
                  // Whether to map submatrices and processors row-wise or column-wise
                  bool isTranspose;
                  
                  
            public : 
            
                  std::pair< int, int > getBlockNumbers(int idSubmatrix) {
                  
                        return blockNumbers[idSubmatrix];
                  
                  }
                  
                  int getTotalSubmatrices() {
                  
                        return totalSubmatrices;
                  
                  }
                  
                  int getColumnsInSubmatrix(int idSubmatrix) {
                  
                        return colSizes[idSubmatrix];
                  
                  }
                  
                  
                  int getColumnIndexStartOfSubmatrix(int idSubmatrix) {
                  
                        return colRanges[idSubmatrix].first;
                  
                  }
                  
                  
                  std::vector< int > getIdSubmatrix(int idProc) {
                  
                        return submatricesPerProc(idProc);
                  
                  }
                  
                  
                  std::vector< int > getRowBreaks() {
                  
                        return rowBreaks;
                  
                  }
                  
                  
                  std::vector< int > getColumnBreaks() {
                  
                        return colBreaks;
                  
                  }
                  
                  
                  std::vector< int > getRowSizes() {
                  
                        return rowSizes;
                  
                  }
                  
                  
                  std::vector< int > getColumnSizes() {
                  
                        return colSizes;
                  
                  }
                  
                  
                  std::pair< int, int > getSubmatrixRowRange(int idSubmatrix) {
                  
                        return rowRanges[idSubmatrix];
                  
                  }
                  
                  
                  std::pair< int, int > getSubmatrixColumnRange(int idSubmatrix) {
                  
                        return colRanges[idSubmatrix];
                  
                  }
                  
                  
                  int getRowsInSubmatrix(int idSubmatrix) {
                  
                        return rowSizes[idSubmatrix];
                  
                  }
                  
                  
                  int getRowIndexStartOfSubmatrix(int idSubmatrix) {
                  
                        return rowRanges[idSubmatrix].first;
                  
                  }
                  
                  
                  int getSubmatrixCountPerRow() { 
                  
                        return submatricesPerRow;
                  
                  }
                  
                  
                  int getSubmatrixCountPerColumn() {
                  
                        return submatricesPerCol;
                  
                  }
                  
                  
                  bool isStoredLocally(int idProc, int idSubmatrix) {
                  
                        return idProc = getIdProcess(idSubmatrix);
                  
                  }
                  
                  
                  void getSubmatrixSize(int idSubmatrix, int & numRows, int & numCols) {
                  
                        numRows = rowSizes[idSubmatrix];
                        numCols = rowSizes[idSubmatrix];
                  
                  }
                  
                  
                  void setTranspose(bool newState) {
                  
                        isTranspose = newState;
                  
                  }
                  
                  
                  bool getTranspose() {
                  
                        return isTranspose;
                        
                  }


                  bool equals(const DistributionBase * other) {
                  
                        if (type != other->type) { return false ; }
                        if (n != other->n) { return false ; }
                        if (m != other->m) { return false ; }
                        if (submatricesPerRow != other->submatricesPerRow) { return false ; }
                        if (submatricesPerCol != other->submatricesPerCol) { return false ; }
                        if (rowBreaks != other->rowBreaks) { return false ; }
                        if (colBreaks != other->colBreaks) { return false ; }
                        if (procPerSubmatrix != other->procPerSubmatrix) { return false ; }
                        if (isTranspose != other->isTranspose) { return false ; }
                        if (offset != other->offset) { return false ; }
                        return true;
                  
                  }
                  
                  
                  int almostEquals(const DistributionBase * other) {
                  
                        if (type != other->type) { return -1 ; }
                        if (n != other->n) { return -1 ; }
                        if (m != other->m) { return -1 ; }
                        if (submatricesPerRow != other->submatricesPerRow) { return -1 ; }
                        if (submatricesPerCol != other->submatricesPerCol) { return -1 ; }
                        if (rowBreaks != other->rowBreaks) { return -1 ; }
                        if (colBreaks != other->colBreaks) { return -1 ; }
                        if (isTranspose != other->isTranspose) { return -1 ; }
                        int result = procPerSubmatrix[0] - other->procPerSubmatrix[0];
                        bool matches = true;
                        for (int submat = 1; submat < totalSubmatrices; ++submat) {
                              if (procPerSubmatrix[submat] - other->procPerSubmatrix[submat] != result) {
                                    matches = false;
                                    break
                              }
                        }
                        if (matches) {
                              return result;
                        } else {
                              return -1;
                        }
                  
                  }

                  
                  int getIdProcess(int idSubmatrix) {
                  
                        return procPerSubmatrix[idSubmatrix];
                  
                  }
                  
                  
                  std::vector< int > getAllIdProcesses() {
                  
                        return procPerSubmatrix;
                        
                  }
                  

                  void getLocalIndices(int rowIndex, int columnIndex, int & idSubmatrix, int & localRowIndex, int & localColumnIndex) { 
                  
                        int rowblock, colblock;
                        if (!(isTranspose)) {
                              for (rowblock = 0; rowblock < submatricesPerCol; ++rowblock) {
                                    if (rowBreaks[rowblock + 1] > rowIndex) {
                                          localRowIndex = rowIndex - rowBreaks[rowblock];
                                          break;
                                    }
                              }
                              for (colblock = 0; colblock < submatricesPerRow; ++colblock) {
                                    if (colBreaks[colblock + 1] > columnIndex) {
                                          localColumnIndex = columnIndex - colBreaks[colblock];
                                          break;
                                    }
                              }
                        } else {                             
                              for (rowblock = 0; rowblock < submatricesPerRow; ++rowblock) {
                                    if (colBreaks[rowblock + 1] > rowIndex) {
                                          localRowIndex = rowIndex - colBreaks[rowblock];
                                          break;
                                    }
                              }
                              for (colblock = 0; colblock < submatricesPerCol; ++colblock) {
                                    if (rowBreaks[colblock + 1] > columnIndex) {
                                          localColumnIndex = columnIndex - rowBreaks[colblock];
                                          break;
                                    }
                              }
                        }
                        idSubmatrix = mapBlocksToSubmatrix(rowblock, colblock);
                                          
                  }


                  int mapBlocksToSubmatrix(int rowblock, int colblock) {
                  
                        if (!(isTranspose)) {
                              return submatrixPerBlock[rowblock][colblock];
                        } else {
                              return submatrixPerBlock[colblock][rowblock];
                        }
                  
                  }
                  
                  
                  void mapSubmatrixToBlocks(int idSubmatrix, int & rowblock, int & colblock) {
                  
                        rowblock = this->blockNumbers[idSubmatrix].first;
                        colblock = this->blockNumbers[idSubmatrix].second;
                        
                  }
                  
                                                      
                  bool transposeId(bool originalIsTranspose, int originalId, int & newId) {
                  
                        if (originalIsTranspose == isTranspose) {
                              newId = originalId;
                              return false;
                        } else {
                              int rowblock, colblock;
                              if (!(isTranspose)) {
                                    rowblock = originalId % submatricesPerCol;
                                    colblock = originalId / submatricesPerCol;
                              } else {
                                    rowblock = originalId / submatricesPerRow;
                                    colblock = originalId % submatricesPerRow;
                              }
                              newId = mapBlocksToSubmatrix(rowblock, colblock);
                              return true;
                        }
                  
                  }
                  
                  
                  std::vector< std::pair< int, int > > cycle(int shift) {
                  
                        std::vector< std::pair< int, int > > result;
                        if (shift % np != 0) {
                              result.resize(totalSubmatrices);
                              for (int submat = 0; submat < totalSubmatrices; ++submat) {
                                    result[submat].first = procPerSubmatrix[submat];
                                    procPerSubmatrix[submat] = (procPerSubmatrix[submat] + shift) % np;
                                    result[submat].second = procPerSubmatrix[submat];
                              }
                              offset = (offset + shift) % np;
                        }
                        return result;
                  
                  }
                  
                  
                  DistributionBase(DistributionBase * original) :
                  
                        type(original->type),            
                        np(original->np),            
                        n(original->n),            
                        m(original->m),            
                        submatricesPerRow(original->submatricesPerRow),            
                        submatricesPerCol(original->submatricesPerCol),            
                        totalSubmatrices(original->totalSubmatrices),            
                        procPerSubmatrix(original->procPerSubmatrix),            
                        submatricesPerProc(original->submatricesPerProc),            
                        offset(original->offset),            
                        rowBreaks(original->rowBreaks),            
                        colBreaks(original->colBreaks),            
                        rowSizes(original->rowSizes),            
                        colSizes(original->colSizes),            
                        rowRanges(original->rowRanges),            
                        colRanges(original->colRanges),            
                        blockNumbers(original->blockNumbers),            
                        submatrixPerBlock(original->submatrixPerBlock),            
                        isTranspose(original->isTranspose)         
                                          
                  {}
                  
                  
                  void populateFromFunctions(int rowblocks, int colblocks, int (*procFunc)(int),
                        int (*rowPartFunc)(int), int (*colPartFunc)(int)) {
                  
                        submatricesPerRow = rowblocks;
                        submatricesPerCol = colblocks;
                        totalSubmatrices = submatricesPerRow * submatricesPerCol;
                        rowBreaks.resize(submatricesPerCol + 1);
                        colBreaks.resize(submatricesPerRow + 1);
                        rowSizes.resize(totalSubmatrices);
                        colSizes.resize(totalSubmatrices);
                        rowRanges.resize(totalSubmatrices);
                        colRanges.resize(totalSubmatrices);
                        blockNumbers.resize(totalSubmatrices);
                        submatrixPerBlock.resize(submatricesPerCol);
                        procPerSubmatrix.resize(totalSubmatrices);
                        submatricesPerProc.resize(this->np);
                        offset = 0;
                        int rowblock, colblock;
                        for (rowblock = 0; rowblock < submatricesPerCol; ++rowblock) {
                              rowBreaks[rowblock] = rowPartFunc(rowblock);
                              submatrixPerBlock[rowblock].resize(submatricesPerRow);
                        }
                        rowBreaks[submatricesPerCol] = this->n;
                        for (colblock = 0; colblock < submatricesPerRow; ++colblock) {
                              colBreaks[colblock] = colPartFunc(colblock);
                        }
                        colBreaks[submatricesPerRow] = this->m;
                        for (int submatrix = 0; submatrix < totalSubmatrices; ++submatrix) {
                              rowblock = submatrix / submatricesPerRow;
                              colblock = submatrix % submatricesPerCol;
                              rowSizes[submatrix] = rowBreaks[rowblock + 1] - rowBreaks[rowblock];
                              colSizes[submatrix] = colBreaks[colblock + 1] - colBreaks[colblock];
                              rowRanges[submatrix] = std::make_pair(rowBreaks[rowblock], rowBreaks[rowblock + 1]);
                              colRanges[submatrix] = std::make_pair(colBreaks[colblock], colBreaks[colblock + 1]);
                              blockNumbers[submatrix] = std::make_pair(rowblock, colblock);
                              submatrixPerBlock[rowblock][colblock] = submatrix;
                              procPerSubmatrix[submatrix] = procFunc(submatrix);
                              submatricesPerProc[procPerSubmatrix[submatrix]].push_back(submatrix);
                        }
                        isTranspose = false;
                        
                  }
 
 
                   void populateFromMaps(std::vector< int > rowcuts, std::vector< int > colcuts, std::vector< int > procs) {
                  
                        submatricesPerRow = rowcuts.size - 1;
                        submatricesPerCol = colcuts.size - 1;
                        totalSubmatrices = submatricesPerRow * submatricesPerCol;
                        rowBreaks.resize(submatricesPerCol + 1);
                        colBreaks.resize(submatricesPerRow + 1);
                        rowSizes.resize(totalSubmatrices);
                        colSizes.resize(totalSubmatrices);
                        rowRanges.resize(totalSubmatrices);
                        colRanges.resize(totalSubmatrices);
                        blockNumbers.resize(totalSubmatrices);
                        submatrixPerBlock.resize(submatricesPerCol);
                        procPerSubmatrix.resize(totalSubmatrices);
                        submatricesPerProc.resize(this->np);
                        offset = 0;
                        int rowblock, colblock;
                        for (rowblock = 0; rowblock < submatricesPerCol; ++rowblock) {
                              rowBreaks[rowblock] = rowcuts[rowblock];
                              submatrixPerBlock[rowblock].resize(submatricesPerRow);
                        }
                        rowBreaks[submatricesPerCol] = this->n;
                        for (colblock = 0; colblock < submatricesPerRow; ++colblock) {
                              colBreaks[colblock] = colcuts[colblock];
                        }
                        colBreaks[submatricesPerRow] = this->m;
                        for (int submatrix = 0; submatrix < totalSubmatrices; ++submatrix) {
                              rowblock = submatrix / submatricesPerRow;
                              colblock = submatrix % submatricesPerCol;
                              rowSizes[submatrix] = rowBreaks[rowblock + 1] - rowBreaks[rowblock];
                              colSizes[submatrix] = colBreaks[colblock + 1] - colBreaks[colblock];
                              rowRanges[submatrix] = std::make_pair(rowBreaks[rowblock], rowBreaks[rowblock + 1]);
                              colRanges[submatrix] = std::make_pair(colBreaks[colblock], colBreaks[colblock + 1]);
                              blockNumbers[submatrix] = std::make_pair(rowblock, colblock);
                              submatrixPerBlock[rowblock][colblock] = submatrix;
                              procPerSubmatrix[submatrix] = procs[submatrix];
                              submatricesPerProc[procPerSubmatrix[submatrix]].push_back(submatrix);
                        }
                        isTranspose = false;
                        
                  }
                  
                  
                  bool compareIdPairs(const std::pair< int, int > & first, const std::pair< int, int > & second) {
                  
                        return (first.second < second.second);
                  
                  }
                 
                  
                  // Both result maps are indexed procId and then submatrixId
                  // sendMap inner first is source procId, second is source submatrixId
                  // recvMap inner first is destination procId, second is destination submatrixId
                  // ASSUMES ONLY DIFFERENCE IS ORDER, NOT PARTITION!!!
                  void producePermutationMap(const DistributionBase * newdist, 
                        std::vector< std::pair< int, int > > & permutationMap); {
                        
                        permutationMap.resize(totalSubmatrices);
                        int curid, destproc, destid;
                        std::vector< int > availableIds;
                        availableIds.resize(totalSubmatrices);
                        std::pair< int, int > cursizes, othersizes;
                        for (int submat = 0; submat < totalSubmatrices; ++submat) {
                              availableIds[submat] = submat;
                        }
                        for (int proc = 0; proc < np; ++proc) {
                              for (int submat = 0; submat < submatricesPerProc[proc].size(); ++submat) {
                                    curid = submatricesPerProc[proc][submat];
                                    cursizes.first = this->rowSizes[curid];
                                    cursizes.second = this->colSizes[curid];
                                    for (int id = 0; id < availableIds.size(); ++id) {
                                          othersizes.first = other->rowSizes[availableIds[id]];
                                          othersizes.second = other->colSizes[availableIds[id]];
                                          if (cursizes == othersizes) {
                                                destid = availableIds.pop(id);
                                                destproc = other->getIdProcess(destid);
                                                permutationMap[curid].first = destid;
                                                permutationMap[curid].second = destproc;
                                                break;
                                          }
                                    }
                              }
                        }
                        
                  }
                  
                  
                  // Tells a single processor how to carve up a single submatrix for a repartitioning
                  void producePartitionSendMaps(const DistributionBase * other, const int idSubmatrix, 
                        std::vector< int > & destProcs, std::vector< int > & destSubmatrices, std::vector< int [4] > & chunks,
                        std::vector< int > & destRbs, std::vector< int > & destCbs) {
                  
                        int rowblock, colblock;
                        this->mapSubmatrixToBlocks(idSubmatrix, rowblock, colblock);
                        int down = this->rowBreaks[rowblock], up = this->rowBreaks[rowblock + 1];
                        int left = this->colBreaks[colblock], right = this->colBreaks[colblock + 1];
                        int u, d, l, r;
                        int chunk [4];
                        for (int rb = 0; rb < this->submatricesPerCol; ++rb) {
                              d = other->rowBreaks[rb];
                              u = other->rowBreaks[rb + 1];
                              for (int cb = 0; cb < this->submatricesPerRow; ++cb) {
                                    l = other->colBreaks[cb];
                                    r = other->colBreaks[cb + 1];
                                    chunk = this->blocksMeet(down, up, left, right, d, u, l, r);
                                    if (chunk[0] != -1) {
                                          destSubmatrices.push_back(other->mapBlocksToSubmatrix(rb, cb));
                                          destProcs.push_back(other->getIdProcess(destSubmatrices[destSubmatrices.size() - 1]));
                                          chunks.push_back(chunk);
                                          destRbs.push_back(rb);
                                          destCbs.push_back(cb);
                                    ;
                              }
                        }                  
                  
                  }
                  
                  // Tells a single processor how to assemble a new submatrix following a repartitioning
                  void producePartitionRecvMaps(const DistributionBase * other, const int idSubmatrix, 
                        std::vector< int > & srcProcs, std::vector< int > & srcSubmatrices, std::vector< int [4] > & chunks) {
                  
                        int rowblock, colblock;
                        other->mapSubmatrixToBlocks(idSubmatrix, rowblock, colblock);
                        int down = other->rowBreaks[rowblock], up = other->rowBreaks[rowblock + 1];
                        int left = other->colBreaks[colblock], right = other->colBreaks[colblock + 1];
                        int u, d, l, r;
                        int chunk [4];
                        for (int rb = 0; rb < other->submatricesPerCol; ++rb) {
                              d = this->rowBreaks[rb];
                              u = this->rowBreaks[rb + 1];
                              for (int cb = 0; cb < other->submatricesPerRow; ++cb) {
                                    l = this->colBreaks[cb];
                                    r = this->colBreaks[cb + 1];
                                    chunk = other->blocksMeet(down, up, left, right, d, u, l, r);
                                    if (chunk[0] != -1) {
                                          destSubmatrices.push_back(this->mapBlocksToSubmatrix(rb, cb));
                                          destProcs.push_back(this->getIdProcess(destSubmatrices[destSubmatrices.size() - 1]));
                                          chunks.push_back(chunk);
                                    ;
                              }
                        }                  
                  
                  }
                  
                  
                  int [4] blocksMeet(int down, int up, int left, int right, int d, int u, int l, int r) {
                  
                        int result [4];
                        bool rowmeet, colmeet;
                        rowmeet = (u > down || d < up);
                        colmeet = (r > left || l < right);
                        if (rowmeet && colmeet) {
                              if (this->n > 1 && this->m > 1) {
                                    result[0] = std::max(down, d);
                                    result[1] = std::min(up, u);
                                    result[2] = std::max(left, l);
                                    result[3] = std::min(right, r);
                              } else {
                                    if (this->isTranspose) {
                                          result[0] = 0;
                                          result[1] = 0;
                                          result[2] = std::max(left, l);
                                          result[3] = std::min(right, r);
                                    } else {
                                          result[0] = std::max(down, d);
                                          result[1] = std::min(up, u);
                                          result[2] = 0;
                                          result[3] = 0;
                                    }
                              }
                        } else {
                              result[0] = -1;
                              result[1] = -1;
                              result[2] = -1;
                              result[3] = -1;
                        }
                        return result;
                  
                  }
                  
                  
                  int [6] getInitData(int idSubmatrix) {
                  
                        int [6] result;
                        result[0] = rowSizes[idSubmatrix];
                        result[1] = colSizes[idSubmatrix];
                        result[2] = rowRanges[idSubmatrix].first;
                        result[3] = colRanges[idSubmatrix].first;
                        result[4] = blockNumbers[idSubmatrix].first;
                        result[5] = blockNumbers[idSubmatrix].second;
                        return result;
                  
                  }
                  
                  
                  bool compareSizePairs(const std::pair< int, int > & first, const std::pair< int, int > & second) {
                  
                        return (first.first < second.first || (first.first == second.first && first.second < second.second));
                  
                  }
                  
                  
                  std::vector< std::pair< int, int > > sortSizes() {
                  
                        std::vector< std::pair< int, int > > result;
                        result.resize(totalSubmatrices);
                        for (int submat = 0; submat < totalSubmatrices; ++submat) {
                              result[submat] = std::make_pair(rowSizes[submat], colSizes[submat]);
                        }
                        std::sort(result, this->compareSizePairs);
                        return result;
                  
                  }
                  
                  
                  bool permutable(DistributionBase * other) {
                        
                        if (n != other->n) { return false ; }
                        if (m != other->m) { return false ; }
                        if (submatricesPerRow != other->submatricesPerRow) { return false ; }
                        if (submatricesPerCol != other->submatricesPerCol) { return false ; }
                        std::vector< std::pair< int, int > > mysortedsizes = this->sortSizes();
                        std::vector< std::pair< int, int > > othersortedsizes = other->sortSizes(); 
                        if (mysortedsizes != othersortedsizes) { return false ; }
                        return true;
                        
                  }
                  
                  
                  bool partitionable(DistributionBase * other) {
                  
                        if (n != other->n) { return false ; }
                        if (m != other->m) { return false ; }
                        return true;                        
                  
                  }
                                  
                  
                  int howDifferent(DistributionBase * other) {
                  
                        int result = 0;
                        if (!(this->equals(other))) {
                              ++result;
                        } else {
                              return result;
                        }
                        if (this->almostEquals(other) == -1) {
                              ++result;
                        } else {
                              return result;
                        }
                        if (!(this->permutable(other))) {
                              ++result;
                        } else {
                              return result;
                        }
                        if (!(this->partitionable(other))) {
                              ++result;
                        }
                        return result;
                  
                  }
                        
                  
                  virtual DistributionBase * clone() = 0;
                  
                  virtual void prepare(int np, int n, int m) = 0;
                                    
                  virtual void populateSelf(int avgr, int avgc) = 0;
                 
                                    
      };
      
      
      class ContiguousDistribution : public DistributionBase {
      
 
            public :                                    
                  
                  
                  ContiguousDistribution(ContiguousDistribution * original) : 
                  
                        DistributionBase(original)
                        
                  {}
                  
                  
                  ContiguousDistribution * clone() {
                  
                        return new ContiguousDistribution(this)(;
                  
                  }
                  
                  
                  void prepare(int np, int n, int m) {
                  
                        this->np = np;
                        this->n = n;
                        this->m = m;
                        type = Contig;
                  
                  }
                  
                  
                  void populateSelf(int avgr, int avgc) {

                        submatricesPerRow = m / avgc;
                        submatricesPerCol = n / avgr;
                        if (m % avgc != 0) { ++submatricesPerRow ; }
                        if (n % avgr != 0) { ++submatricesPerCol ; }
                        totalSubmatrices = submatricesPerRow * submatricesPerCol;
                        rowBreaks.resize(submatricesPerCol + 1);
                        colBreaks.resize(submatricesPerRow + 1);
                        rowSizes.resize(totalSubmatrices);
                        colSizes.resize(totalSubmatrices);
                        rowRanges.resize(totalSubmatrices);
                        colRanges.resize(totalSubmatrices);
                        blockNumbers.resize(totalSubmatrices);
                        submatrixPerBlock.resize(submatricesPerCol);
                        procPerSubmatrix.resize(totalSubmatrices);
                        submatricesPerProc.resize(this->np);
                        offset = 0;
                        int rowblock, colblock;
                        for (rowblock = 0; rowblock < submatricesPerCol; ++rowblock) {
                              rowBreaks[rowblock] = rowblock * avgr;
                              submatrixPerBlock[rowblock].resize(submatricesPerRow);
                        }
                        rowBreaks[submatricesPerCol] = this->n;
                        for (colblock = 0; colblock < submatricesPerRow; ++colblock) {
                              colBreaks[colblock] = colblock * avgc;
                        }
                        colBreaks[submatricesPerRow] = this->m;
                        for (int submatrix = 0; submatrix < totalSubmatrices; ++submatrix) {
                              rowblock = submatrix / submatricesPerRow;
                              colblock = submatrix % submatricesPerCol;
                              rowSizes[submatrix] = rowBreaks[rowblock + 1] - rowBreaks[rowblock];
                              colSizes[submatrix] = colBreaks[colblock + 1] - colBreaks[colblock];
                              rowRanges[submatrix] = std::make_pair(rowBreaks[rowblock], rowBreaks[rowblock + 1]);
                              colRanges[submatrix] = std::make_pair(colBreaks[colblock], colBreaks[colblock + 1]);
                              blockNumbers[submatrix] = std::make_pair(rowblock, colblock);
                              submatrixPerBlock[rowblock][colblock] = submatrix;
                              procPerSubmatrix[submatrix] = ((rowblock * submatricesPerRow) + colblock) / this->np;
                              submatricesPerProc[procPerSubmatrix[submatrix]].push_back(submatrix);
                        }
                        isTranspose = false;
                        
                  }

      };


      class CyclicDistribution : public DistributionBase {
      
 
            public :
                   
                  CyclicDistribution(CyclicDistribution * original) : 
                  
                        DistributionBase(original)
                        
                  {}
                  
                  
                  CyclicDistribution * clone() {
                  
                        return new CyclicDistribution(this)(;
                  
                  }
                  
                  void prepare(int np, int n, int m) {
                  
                        this->np = np;
                        this->n = n;
                        this->m = m;
                        type = Cyclic;
                  
                  }
                  
                  
                  void populateSelf(int avgr, int avgc) {

                        submatricesPerRow = m / avgc;
                        submatricesPerCol = n / avgr;
                        if (m % avgc != 0) { ++submatricesPerRow ; }
                        if (n % avgr != 0) { ++submatricesPerCol ; }
                        totalSubmatrices = submatricesPerRow * submatricesPerCol;
                        rowBreaks.resize(submatricesPerCol + 1);
                        colBreaks.resize(submatricesPerRow + 1);
                        rowSizes.resize(totalSubmatrices);
                        colSizes.resize(totalSubmatrices);
                        rowRanges.resize(totalSubmatrices);
                        colRanges.resize(totalSubmatrices);
                        blockNumbers.resize(totalSubmatrices);
                        submatrixPerBlock.resize(submatricesPerCol);
                        procPerSubmatrix.resize(totalSubmatrices);
                        submatricesPerProc.resize(this->np);
                        offset = 0;
                        int rowblock, colblock;
                        for (rowblock = 0; rowblock < submatricesPerCol; ++rowblock) {
                              rowBreaks[rowblock] = rowblock * avgr;
                              submatrixPerBlock[rowblock].resize(submatricesPerRow);
                        }
                        rowBreaks[submatricesPerCol] = this->n;
                        for (colblock = 0; colblock < submatricesPerRow; ++colblock) {
                              colBreaks[colblock] = colblock * avgc;
                        }
                        colBreaks[submatricesPerRow] = this->m;
                        for (int submatrix = 0; submatrix < totalSubmatrices; ++submatrix) {
                              rowblock = submatrix / submatricesPerRow;
                              colblock = submatrix % submatricesPerCol;
                              rowSizes[submatrix] = rowBreaks[rowblock + 1] - rowBreaks[rowblock];
                              colSizes[submatrix] = colBreaks[colblock + 1] - colBreaks[colblock];
                              rowRanges[submatrix] = std::make_pair(rowBreaks[rowblock], rowBreaks[rowblock + 1]);
                              colRanges[submatrix] = std::make_pair(colBreaks[colblock], colBreaks[colblock + 1]);
                              blockNumbers[submatrix] = std::make_pair(rowblock, colblock);
                              submatrixPerBlock[rowblock][colblock] = submatrix;
                              procPerSubmatrix[submatrix] = ((rowblock * submatricesPerRow) + colblock) / this->np;
                              submatricesPerProc[procPerSubmatrix[submatrix]].push_back(submatrix);
                        }
                        isTranspose = false;
                        
                  }
                        
      };


      class FreeDistribution : public DistributionBase {
      
 
            public :
                   
                  FreeDistribution(FreeDistribution * original) : 
                  
                        DistributionBase(original)
                        
                  {}
                  
                  
                  FreeDistribution * clone() {
                  
                        return new FreeDistribution(this)(;
                  
                  }
                  
                  void prepare(int np, int n, int m) {
                  
                        this->np = np;
                        this->n = n;
                        this->m = m;
                        type = Free;
                  
                  }
                  
                  
                  // Cannot populate self, requires maps or functions
                  void populateSelf(int avgr, int avgc) {

                        return;
                        
                  }

      };


}
